﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace TimbreRecognition.FeatureExtraction
{
   /*
      Please feel free to use/modify this class. 
      If you give me credit by keeping this information or
      by sending me an email before using it or by reporting bugs , i will be happy.
      Email : gtiwari333@gmail.com,
      Blog : http://ganeshtiwaridotcomdotnp.blogspot.com/ 
   */

    public class MelCepstrum {

        private const double LowerFilterFreq = 0;

        private double UpperFilterFreq = 22050;

	    private int numMelFilters = 30;// 24 is standard?
	    private int numCepstra;// number of mfcc coeffs
	    private double preEmphasisAlpha = 0.95;

	    private double samplingRate;
	    private double[] bin;
	    private int samplePerFrame;
	    // /////
	    private FFT fft;
	    private DCT dct;

        public MelCepstrum(/*int samplePerFrame,*/ int samplingRate, int numCepstra)
        {
		    //this.samplePerFrame = samplePerFrame;
		    this.samplingRate = samplingRate;
		    this.numCepstra = numCepstra;
		    //upperFilterFreq = samplingRate / 2.0;
		    fft = new FFT();
		    dct = new DCT(this.numCepstra, numMelFilters);
	    }

	    public double[] doMFCC(float[] framedSignal) {
		    // Magnitude Spectrum

            samplePerFrame = framedSignal.Length;

		    bin = magnitudeSpectrum(framedSignal);
		    framedSignal = preEmphasis(framedSignal);
		    /*
		     * cbin=frequencies of the channels in terms of FFT bin indices (cbin[i]
		     * for the i -th channel)
		     */

		    // prepare filter for for melFilter
		    int[] cbin = fftBinIndices();// same for all
		    // process Mel Filterbank
		    double[] fbank = melFilter(bin, cbin);
		    // magnitudeSpectrum and bin filter indices

		    // System.out.println("after mel filter");
		    // ArrayWriter.printDoubleArrayToConole(fbank);

		    // Non-linear transformation
		    double[] f = nonLinearTransformation(fbank);
		    // System.out.println("after N L T");
		    // ArrayWriter.printDoubleArrayToConole(f);
		
		    // Cepstral coefficients, by DCT
		    double[] cepc = dct.performDCT(f);
		    // System.out.println("after DCT");
		    // ArrayWriter.printDoubleArrayToConole(cepc);
           
            double[] main = new double[cepc.Length];

	        Array.Copy(cepc, main, main.Length);

	        return main;
	    }

	    private double[] magnitudeSpectrum(float[] frame) {
		    double[] magSpectrum = new double[frame.Length];
		    // calculate FFT for current frame
		    fft.computeFFT(frame);
		    // System.err.println("FFT SUCCEED");
		    // calculate magnitude spectrum
		    for (int k = 0; k < frame.Length; k++) {
			    magSpectrum[k] = Math.Sqrt(fft.real[k] * fft.real[k] + fft.imag[k] * fft.imag[k]);
		    }
		    return magSpectrum;
	    }

	    /**
	     * emphasize high freq signal
	     * 
	     * @param inputSignal
	     * @return
	     */
	    private float[] preEmphasis(float[] inputSignal) {
		    // System.err.println(" inside pre Emphasis");
		    float[] outputSignal = new float[inputSignal.Length];
		    // apply pre-emphasis to each sample
		    for (int n = 1; n < inputSignal.Length; n++) {
			    outputSignal[n] = (float) (inputSignal[n] - preEmphasisAlpha * inputSignal[n - 1]);
		    }
		    return outputSignal;
	    }

	    private int[] fftBinIndices() {
		    int[] cbin = new int[numMelFilters + 2];
		    cbin[0] = (int) Math.Round(LowerFilterFreq / samplingRate * samplePerFrame);// cbin0
            cbin[cbin.Length - 1] = (int)Math.Round(UpperFilterFreq / samplingRate * samplePerFrame);// cbin24S
		    //cbin[cbin.Length - 1] = (samplePerFrame / 2);// cbin24
		    for (int i = 1; i <= numMelFilters; i++) {// from cbin1 to cbin23
			    double fc = centerFreq(i);// center freq for i th filter
			    cbin[i] = (int) Math.Round(fc / samplingRate * samplePerFrame);
		    }
		    return cbin;
	    }

	    /**
	     * performs mel filter operation
	     * 
	     * @param bin
	     *            magnitude spectrum (| |)^2 of fft
	     * @param cbin
	     *            mel filter coeffs
	     * @return mel filtered coeffs--> filter bank coefficients.
	     */
	    private double[] melFilter(double[] bin, int[] cbin) {
		    double[] temp = new double[numMelFilters + 2];
		    for (int k = 1; k <= numMelFilters; k++) {
			    double num1 = 0.0, num2 = 0.0;
			    for (int i = cbin[k - 1]; i <= cbin[k]; i++) {
				    // System.out.println("Inside filter loop");
				    num1 += ((i - cbin[k - 1] + 1) / (cbin[k] - cbin[k - 1] + 1)) * bin[i];
			    }

			    for (int i = cbin[k] + 1; i <= cbin[k + 1]; i++) {
				    // System.out.println("Inside filter loop 222222");
				    num2 += (1 - ((i - cbin[k]) / (cbin[k + 1] - cbin[k] + 1))) * bin[i];
			    }

			    temp[k] = num1 + num2;
		    }
		    double[] fbank = new double[numMelFilters];
		    for (int i = 0; i < numMelFilters; i++) {
			    fbank[i] = temp[i + 1];
			    // System.out.println(fbank[i]);
		    }
		    return fbank;
	    }

	    /**
	     * performs nonlinear transformation
	     * 
	     * @param fbank
	     * @return f log of filter bac
	     */
	    private double[] nonLinearTransformation(double[] fbank) {
		    double[] f = new double[fbank.Length];
		    const double FLOOR = -50;
		    for (int i = 0; i < fbank.Length; i++) {
			    f[i] = Math.Log(fbank[i]);
			    // check if ln() returns a value less than the floor
			    if (f[i] < FLOOR) {
				    f[i] = FLOOR;
			    }
		    }
		    return f;
	    }

	    private double centerFreq(int i) {
		    double melFLow, melFHigh;
		    melFLow = freqToMel(LowerFilterFreq);
		    melFHigh = freqToMel(UpperFilterFreq);
		    double temp = melFLow + ((melFHigh - melFLow) / (numMelFilters + 1)) * i;
		    return inverseMel(temp);
	    }

	    private double inverseMel(double x) {
		    double temp = Math.Pow(10, x / 2595) - 1;
		    return 700 * (temp);
	    }

	    protected double freqToMel(double freq) {
		    return 2595 * log10(1 + freq / 700);
	    }

	    private double log10(double value) {
		    return Math.Log(value) / Math.Log(10);
	    }
    }
}
